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Abstract 

The new method of solution of the Bethe-Salpeter equation for quark-antiquark pseudoscalar 
bound state is proposed. With the help of integral representation the results are directly obtained 
in Minkowski space. Dressing of Greens functions is naturally considered providing thus the correct 
inclusion of the running coupling constant and the quark propagators as well as. The first numerical 
results are presented for a simplified ladder approximation. 

PACS numbers: 11. 10. St, 11.15.Tk 



1 



I. INTRODUCTION 



Among various approaches used in meson physics the formalism of Bethe-Salpeter and 
Dyson- Schwinger equations (DSE's) plays the traditional and indispensable role. Bethe- 
Salpeter equation (BSE) provides quantum field theoretical starting point to describe 
hadrons as relativistic bound states of quarks (antiquarks). DSE's and BSE's framework 
has been widely used in order to obtain the nonperturbative information about the spec- 
tra and decays of the whole lowest pseudoscalar nonet with emphasizing attention to the 
QCD pseudo-Goldstone bosons - the pions [l|. With certain success the formalism offers 
a satisfactory window into the 'next scale' meson sector including also excited and vector, 
scalar mesons 0]. Within formalism, the mesons electromagnetic form factors have been 
calculated at spacelike regime of momenta [3|. 

When dealing with bound states composed from the light quarks then there is unavoidable 
necessity of the full use of BSE covariant framework. The nonperturbative knowledge of the 
Greens function composing the kernel of BSE is required. Very often, the problem is solved 
in the Euclidean space where it is more tractable since there is no singularities in Green 
functions there. The physical amplitudes can be then obtained by continuation to the 
Minkowski space. Note, the extraction of the mass spectra is already a complicated task 
[H, while there is a practical absence of an analytical continuation of Euclidean space form 
factors. 

When dealing with the heavy quarkonia or with the mesons composed from the mix 
like B c meson (the state measured in the Fermilab Tevatron by CDF Collaboration) some 
simplifying approximations are possible. Different approaches have been developed to reduce 
computational complexity of full 4-dim BSE. The so called instantaneous approximation [B| 
or quasi-potential approximation 0] can reduce the 4-dim BSE to the 3-dim equation in 
a Lorentz covariant manner. In practise, such 3-dim equations are more tractable since 
the solution is less involved, especially when one tricky exploits a considerable freedom in 
performing three dimensional reduction. Note, contrary to ladder BSE, these equations also 
reduce to the Schrodinger equation of non-relativistic Heavy Meson Effective Theory and 
nonrelativistic QCD 0]. However, the interaction kernels of the reducted equations often 
represents the input based on economical model phenomenology and the connection to the 
underlying theory (QCD) is less clear (if not abandoned from the beginning). 

In this paper, we extend the method of solution of the full 4-dim BSE, originally devel- 
oped for the pure scalar theories [s| , [l(| , [H| , to the theories with non-trivial spin degrees of 
freedom. Within a certain assumptions on the functional form of Greens functions we de- 
velop the method of direct Minkowski space BSE solution in its original manifestly Lorentz 
covariant (four-dimensional) manner. In order to make this paper as self-contained as pos- 
sible, we will supply some basic facts concerning the BSE approach to relativistic meson 
bound states bellow here. 

The crucial step to derive the homogeneous BSE for bound state is the assumptions that 
the bound state reflects itself in a pole in the four-point Green's function for an on-shell 
momentum P, P 2 = M 2 



GW(p,f/, P)=E , ^ P os)Wi Pes) (L1) 

' y (2n)t2E Pj (P°-E p . +ie) S 1 ; 



where E Pj = yjp + Mj, Mj is a positive mass of the bound state characterized by BS wave 
function ipj within a set of quantum numbers j . Then the BSE can be conventionally 
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written in momentum space like: 



5r 1 (p + , P)V>(p, P)5' 2 - 1 (p_, P) = -i J 70j^(p, fc, P)V>(p, ^) , (1-2) 

p+ = p + aP , 

P- = p — (1 — a)P , 



or equivalently written in the term of BS vertex function V: 

T(p, P) = -if -0r 4 V(p, k, P)Sx(k+, P)r(p, P)S 2 (k^ P) , (1.3) 

where we suppress all the Dirac, flavor and Lorentz indices and a G (0, 1). The function V 
represents two body irreducible interaction kernel and Si, i = 1,2 are the dressed propagators 
of the constituents. The free propagators reads 

s?(p) = • . (i-4) 

pZ _ m z _|_ le 



Concerning the solution of BSE fll.3p for pseudoscalar meson, it has the generic form [12 



F{q, P) = 75 [rA + T B q.P 4 + T c P + T D i P + T E P ft (1.5) 

where i — A, B, C, D, E are the scalar functions of their arguments P, q. If the bound 
state has well defined charge parity, say C = 1 then they are even function of q.P and 
Td = — r^- 

As was already discussed in the paper [13j, the dominant contribution to the BSE cal- 
culation for pseudoscalar meson comes namely from the first term in Eq. (11.51) . Within 
15 percentage accuracy it is already true for the light mesons (jr,K,rj) while in the case of 
heavy meson pseudoscalar ground state (r] c , rjb) the contributions due to the other tensor 
components in Eq. (11.51) are even more negligible. Hence, in this stage of our Minkowski 
calculation we also approximate our solution by taking T = j^Ta- 

The interaction kernel is approximated by the dressed gluon propagator with the inter- 
action gluon-quark-antiquark vertices taken in their bare forms, thus we can write 

V(p, q, P) = g 2 (K)D^(p - q, k)^ (g) Y , (1-6) 

where the full gluon propagator is renormalized at the scale k and the effective strong running 
coupling are related through the following equations: 

9 \k)D, v {1,k) =a.(/,«)^-^(«) 7 rrV, (1-7) 



P + ie ^ v 'l 4 + ie 



a s (q, k) 



Ii{q\K) 



From the class of £ linear covariant gauges the Landau gauge £ = is employed thorough 
the presented paper. 
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In the next Section we will derive the solution for the dressed ladder approximation of 
BSE, i.e. the the all propagators are considered as a dressed ones and no crossed diagrams 
are taken into the calculation. The BSE for quark- ant iquark states has been many times 
considered in Euclidean space even beyond ladder approximation. Notably, the importance 
of dressing the proper vertices in the light quark sector has been already stressed in 14 



so 



our approximations is known to be certainly limited. Going beyond rainbow (7^ approxima- 
tion) is straightforward but rather involved (for comparison, see the SDE Minkowski study 
published in 15 1 and in [lfjj, wherein the later study includes the minimal gauge covariant 
vertex instead of bare one). In this paper we prefer to explain the computational method 
instead of making a BSE study with the most sophisticated kernel known in the literature. 

The layout of the papers is as follows: In the next section we describe the method of the 
solution. As an exhibition of the method, the numerical results are presented in the Section 
3. We conclude in Section 4. 



II. INTEGRAL REPRESENTATION AND THE SOLUTION OF BSE 

In this section we describe our method of the BSE solution in Minkowski space. It 
basically assumes that the various Greens functions that appears in the interaction kernel 
of BSE can be written as a weighted integral over the various spectral functions (the real 
distribution) p. 

More explicitly stated, the full quark and gluon (in Landau gauge) propagators are as- 
sumed to satisfy the standard Lehmann representation, which reads 

s{i)= r d ^y +p ^\ (2.1) 

Jo V — uj + ie 



0^(0= f ^/^L O), ( 2 - 2 ) 



I 2 — to + ie 

The representations [2TT| 12.21 are strictly related with the singularity structure and analyt- 
ical properties of functions used for. Until now, within a certain limitations, these integral 
representations have been used for nonperturbative evaluation of Greens functions in vari- 
ous models [13], notin g he re, the true analytical structure of QCD Greens functions is not 
answered (see also [3, 20]) with reliability. 



Alternatively, the integral representation can be more complicated and generalized to the 
product of such weighted integrals. In fact, our treatment shall counts more general integral 
representation for quark propagator like the one recently considered in the paper 22|. Also 



the following alternative but notable possibility for gluon propagator, which read now: 



will be automatically included in our treatment. Note, this is effectively the gluon form factor 
rather then gluon propagator itself which satisfies the appropriate integral representation 



in this case. Such expression, motivated by analyticized running coupling method 18], has 



been recently considered in the studies of dynamical chiral symmetry breaking in various 



approximations [21|, [22], [23 



Further, here we generalize the idea of Perturbation Theory Integral Representation 
(PTIR) [s| to our case. The PTIR represents unique integral representaion for n-point 
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Green function defined by n-leg Feynman integral. The inclusion of dressed propagator rep- 
resenting by the internal lines of (skeleton) Feynman diagram is selfconsistently achieved by 
their own IR (see fl2.2l srforgluon)). 

The generalized PTIR formula for n-point function in theory involving fields with ar- 
bitrary spin is exactly the same as in the original scalar theory considered in [s| but the 
spectral function receives nontrivial tensor structure. Let us denote p(a, Xi) such generalized 
weight function. Then, it can be clearly decomposed into the following sum 

p(a,Xi) scaLt h. -> ^pjia^^Vj (2.4) 
j 

where a,xi represent the set of spectral variables and where j runs over the all possible 
independent combinations of Lorentz tensors and Dirac matrices Pj. The function pj(a,Xi) 
represents just the PTIR weight function of j— form factor (the scalar function by its 
definition) since it can clearly written as a suitable scalar Feynman integral. Leaving the 
question of (renormalization) scheme dependence, we refer the reader to Nakanishi book [s| 
for a detailed derivation of PTIR. 

Let us apply our idea to the pseudoscalar bound state vertex function. Recognizing that 
the singularity structure (given by denominators) of the rhs. of the BSE is the same as in 
the scalar models studied in [Tl[ the appropriate IR for pseudoscalar bound state vertex 
function r^(g, P) has to read 

T A (q, P) = r duo f 1 dz P^^A (2.5) 
Jo J-i [F(u, z; P, q)] 

where we have introduced useful abbreviation for the denominator of the IR (12.5ft 

F(u, z;P,q) = co- (q 2 + q.Pz + P 2 / 4 ) - »"e > (2-6) 

where N is an integer free parameter. The same IR could apply for the other components 
of the full function V which however will not be considered furthermore. 

Substituting IR's (I2.5p . (l2.2p . (l2.ip into the rhs. of BSE 11.31 one can analytically integrate 
over the loop momenta. Assuming Theorem of uniqueness [8| we should arrive to the same 
IR ( 12.51) because of the rhs. of BSE 11.31 The derivation is presented in the Appendix A for 
the cases N = 1,2. 

In the other words, we have converted the momentum BSE (with singular kernel) to the 
homogeneous two dimensional integral equation for the real weight function p^(u,z) 

pjf 1 (oo, z) = duo f 1 dzV [N] (Co, z; oo, z)p [ * ] (oo, z) (2.7) 
Jo J-l 

where the kernel V^(oo, z; to, z) is the regular multivariable function. 

The kernel also automatically supports the domain Q where the function p^{u,z) 
is nontrivial. This domain is always smaller then the infinity strip 0, oo x — 1, 1 as is explic- 
itly assumed by the integral boundaries of oo, z integrals. For instance, for simplest kernel 
parametrized by the free gluon propagator and with constituent quarks of mass m, then for 
the flavor singlet meson we have p^\uo, z) ^ only for oo > m 2 . 

In our approach, to solve the momentum BSE in Minkowski space is equivalent to the 
finding of the real solution of the real integral equation (12. 7p . No special choice of the frame 
or has been used. If one needs the resulting vertex function can be obtained by the numerical 
integration over the p^ in an arbitrary reference frame. 
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III. NUMERICAL RESULTS 



In this section we discuss a numerical solutions of the BSE with a various interacting 
kernels. For this purpose we will vary the coupling strength and the effective gluon mass 
m g as well as. We mainly concern on the range of binding energy that coincides with the 
heavy quarkonia - the object of future interest of us. Also we take a discrete mass m g 
such that it runs from zero to the value of constituent mass of quarks, these values are 
expected to be relevant for the case of the real gluon propagator (when m g is replaced 
by the continuous spectral variable u (12.2)) ) . Thus, in each case, the corresponding gluon 



density is p g (c) = N g 5(c 



which specify the kernel of BSE to be (in Landau gauge): 



V(q-p) = g 2 



(<i-p)n(q-p)v 
(<?-p) 2 



(q-pf 



m. 



IE 



-7 ® T 



(3.1) 



where the prefactor (including the trace of color matrices) is simply absorbed into the 
coupling constant. For our actual calculation we have used the bare constituent propagators 
Si(pi) with heavy quark mass M (see the Appendix A for this approximation). 

Firstly, we followed the standard procedure: after fixing the bound state mass (P 2 ) we 
look for the solution by iterating the BSE for spectral function with fixed coupling constant 
a = g 2 /(4n). Very similar to the scalar case the choice N = 2 of the power of F in 
the IR for the bound state vertex function is preferable one. It reasonably make provision 
for a numerical errors from one side and the computational obstacles for high N from the 
other side, noting here that using N = 1 is rather unhelpful (comparing with massive 
Wick-Cutkosky model) since we did not find the stable solution for a wide class of input 
parameters g, m g for that choice. On the other hand, using the value N = 2, we have found 
the results for the all possible interaction kernels considered here. It also includes the cases 
with vanishing m g , which means the numerical troubles originally presented in the scalar 
models are fully overcame here. The details of our numerical treatment are presented 
in Appendix B. 

As it is more usually in the non-relativistic case, we fix the coupling constant a = g 2 / (Air) 
and then look for the bound state mass spectrum. We have found the same result in the both 
cases, whether we fixed P or a as a first, noting that in the later case the whole integration 
in the kernel K needs to be performed at each iteration step, which makes the problem more 
computer time consuming. 

The obtained solutions for a's for given fractional binding r/ = \/p~ 2 ~/2m = 0.95 and 
various masses m g are displayed in the Table I. When we have fixed the gluon mass as 
m g = 0.5 and varied the fractional binding binding rj we have obtained the spectra as they 
are displayed in the Table II. 



m g /m q 


1(T 3 


0.01 


0.1 


0.5 


a : 


0.666 


0.669 


0.745 


1.029 



TABLE 1. Coupling constant a s = g 2 /(4ir) for several selections of m g /m q 



(m q is the 



constituent quark mass) with given fraction of binding rj = yP 2 /2m = 0.95. 



V ■ 


0.8 


0.9 


0.95 


0.99 


a 


1.20 


1.12 


1.03 


0.816 



6 



0.8 
0.7 



0.6 




FIG. 1: The rescaled weight function r = p jff for the following parameters of the model: 
i] = 0.95, m g = O.OOlMq, a s = 0.666, the small mass m g approximates the one gluon exchange 
interaction kernel. 

TABLE 2. Coupling a s = g 2 / (4ir) as a function of fraction of binding r] = y/P^/2m for 
exchanged massive gluon with m g = 0.5m q , m q is the constituent quark mass. 

For the illustration the weight function jffl is displayed in Fig. ([T]). 

IV. SUMMARY AND CONCLUSIONS 

The main result of this paper is development of technical framework for the solution of 
Bethe-Salpeter equation in Minkowski space. To obtain the spectrum, no preferred reference 
frame is needed and the wave function can be obtained in an arbitrary frame (without 
numerical boosting) by simple integration of the weight function. 

The treatment is based on the usage Integral representation for Green's function of given 
theory, including the bound state vertices themselves. The method has been explained and 
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checked numerically on the samples of fermion-antifermion pseudoscalar bound state. It was 
shown that the momentum space BSE can be converted to the real equation for the real 
weight functions p which is easily solvable numerically. The main motivation of the author 
was to develop the practical tool that respect the self-consistence of Schwinger-Dyson and 
Bethe-Salpeter equations. Generalizing this study to the other mesons (vectors, scalars), 
considering more general flavor (izospin) structure with the simultaneous improvement of 
the approximations (correctly dressed gluon propagator, dressed vertices, etc.) will be a 
basic steps towards the fully Lorentz covariant description of plethora of transitions and 
form factors in timelike region of fourmomenta (Minkowski space). 



APPENDIX A: KERNEL FUNCTIONS 



Writing explicitly Dirac indices the BSE for quark-antiquark bound state reads 

r(g, P) wp = i J -0^S(q + P/2) /37 r(g, P)^S(q - P/2) Y ^V w ^ p (q, k; P) (Al) 

where the Lorentz indices of the vertex function has not been specified. 

In our approximation the integral representation for pseudoscalar bound state vertex is 

r (! ,P)^ T f/>/^,/!M ? ( A2) 

Jo J -i [F(uj,z;P,q)} 
and the generalized kernel for ladder approximation of the BSE reads 

, 2 r a~ 9 » v + (z-p) 2 



vw^Vo dc (g -p )a -V+fc *< c >frft (A3) 

where a, /3, 7, 8(p,, v) indices show up the appropriate Dirac (Lorentz) structure. In the 
addition we use the IR for the all functions entering the BSE, this includes (IA2IIA3j) and the 
ones for propagators S(q±) (12.11) . 

For purpose of brevity we will use the following abbreviation for the prefactor: 

/*oo /*oo /*oo 

du dz dc da dbp [N \uo, z)g 2 p g (c) . (A4) 



-1 

Within this convention the BSE can be written like 



/ 



dujdz- 



p N (u,~z) 



[F(£j,z;p,P)} 



N 



d 4 q p v (a)p v (b)(q 2 - P 2 /4) - p s (a)p s (b) 
(2^) 4 



[F(oo,z;q,P)f D X D 2 D 



where the trace was taken over the Dirac indices (after multiplying by 75) and 



D 3 



+ P/2) 2 - a + ie. 
-P/2f-b + ie, 



(q-pf 



c + ie . 



(A5) 



(A6) 



In what follows, we will transform the rhs. of BSE flA5j) into the integral representation, ie. 
lhs. of BSE flA5l) . 



S 



As a first we use the following algebraic identities: 



_u- q.Pz - P 2 /4 
F(u,z;q,P) ~ F(uo,z;q,P) 

q.P 1/1 1 b-a s 
D X D 2 ~ 2 \D X ~ D~ 2 + D 1 D 2/ 

which gives us for the rhs. of (1A5I) : 



(A7) 
(A8) 



d A q \ -fp„(a)p„(6) 

N 



1 



+ 



sJ (2tt) 4 l[F(w, Z ; 9 ,P)f D 3 VA 7J 

- p2 /2 + °dl z ) pv ( a ) pv ( b ) _ p 5 (q)p s (&) 

[F(w, z; g, P)] N D 1 D 2 D 3 



p v (a)p v (b) 



[F(uj,z;q,P)} N - 1 D 1 D 2 D 3 



(A9) 



In addition we will use Feynman parametrization, starting with the first term in the Rel. 
(IA9I) we can get: 



[F^^P)]" VA -D 



r(iv) 



g 2 + q.P(zx - (1 - x)) + P 2 /A - ux — a(l - x) 
q 2 + q.P(zx + (1 - x)) + P 2 /4 ~ux- 6(1 - x) 



-N-l 



-N-l 



(A10) 



Substituting (1A10I) back to the relation (|A9j) . using the Feynman variable ?/ for matching 
with the scalar propagator D 3 and integrating over the four-momentum q we get for the first 
line in (IA9I) the following expression: 



jv-i f Pv{a)pv{b) f 1 , f 1 





sJV-l 












P 2 r 






T 


y- 




p2 






4 


y 



s 2(4vr) 2 



/ dx [ dyy N x l - N z 
Jo Jo 



-i -AT 



2/ - y 2 (x(l + z) - l) 2 + g 2 (l - y)y + Ppy(l - y){x(l + z) - 1) - P(a) 



y - y 2 (x(^ - 1) + l) 2 + g 2 (l - y)y + Ppy(l - y)(x(z - 1) + 1) - U(b) 



-N ' 



(Ah) 



where we have labeled 

17(a) = [ux + a(l - x)] y + c(l - y) . (A12) 

Taking the substitution x — > z such that z = x(l + z) — lin the firts line inside the braces 
we can write (including prefactors): 



Is 



(4tt) 



7o ^/-Ml 



1 + 5 



~ \ AT-1 



p v {a)p v (b) 



l + z 



2(1 -y) N F(£l,z;p,P) 



N ' 



(A13) 



where F is defined by ( 12. 6ft and 



fg + qf-^ + c(l-y)-^ 2 (l-z 2 ) 
2/(1-2/) 



(A14) 



9 



Introducing the identity 

/■oo _, 

dw8{Cb - O) (A15) 



into the Rel. (IA14I) . changing the integral ordering and integrating over the Feynman vari- 
able y we can obtain the desired expression: 



du I dz ^ ' ' 



i [F(a,,i;p,P)r 



Xl(a; ' z) -"isw i+* (i-^-r- 1 ^ ' (A16) 

where y a ± are the roots of the quadratic equation 

y 2 A + y B a + c = (A17) 
with the functions A, B a , D a defined like 

A = Co - S; B a = (lo- a)T + + a - c - uj; D a = B 2 a - AAc 

P 2 1 + 7 

S = (l-z 2 ) — ; T± = — . (A18) 
V ; 4 ' l±z V 7 

Repeating the similar procedure in the second term in rel. ( 1A11I) we can write for this: 
v (u ~z) - ( ^ Z o (a)o (b)9(S z) V ~ Mg) 

y 6± = B 6 = (a; - &)T_ + 6 - c - u; ;D b = B 2 - AAc (A19) 

where we used the label %2 instead of xi- 

In what follows we transform the last term of the second line of Rel. (1A9j) to the desired 
form of integral representation (12.51) . For this purpose we basily follow the derivation already 
presented in since the procedure is exactly the same. However, due to notational 
difference we present all calculational details bellow. 

Let us denote the relevant integral as: 



i ,nJ\A~ „ „ ^iv^T (A20) 



d 4 q 1 
(2n^ Dl D 2 D 3 [F{u,z;q,P)] 
Using the Feynman parameterization technique we first write 



D ^ = 2 / 

-1 



If drj 



[M 2 — f(q, P, rj) — it 



2 ■ 



M 2^a + b + ^_b v . f ^p^ ) = q 2 + r]q .p + P 2 /A . (A21) 
Next the denominator of IR for bound state vertex is added: 



l l 



D 1 D 2 r(JV + l) f, t ^ (1 - t)t 



N-2 



[F(u,z;q,P)] N ^ 2T(N - 1) 7 'J [R — f(q, P, z) — ie] N+1 
R = ujt+(l-t)M 2 , (A22) 
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where z = tz + (1 — t)r]. Now, we match with the function D3 and integrate over the four 
momentum q. Thus we get 



/ = % 



d A q DiD 2 Dz 
(2n)*iF(u;,z- q ,P)] N -i 

r(iV + 2) Jd V Jdt(l-t)t N - 2 Jdxx N I, 



(A23) 



—1 



2V{N -I) 



1 -i 



-1 



WT 4 



-Q 1 + Q ■ Q - (1 - x)p 2 - -P 2 + (1 - x)c + xR-ie 



-(N+2) 



r(N) 



1 



1 



n(t) 



(4vr) 2 r(iV + 2) 2^(1 - z)* - f(q, P, z) - ie] N 

R c x 
+ - 



1 — x x (1 — x] 



S. 



(A24) 



where Q = (1 — x)p — xzP/2. The function S" is defined by (1A18I) and note z lies in the 
interval < —1, +1 >, < S < (y/a + \/b) 2 /4. Interchanging the integrals over rj and t such 
that : 



! 1 1 

drj dt = dz 



-1 



T + dt 



l-t 



Q(z-z) + 



T ~ dt 



l-t 



e{z-z) 



- ± 



1 ±5 

r±~z 



and z = tz + {1 — t)rj , 



we finally obtain 



iV-1 
2(4vr) 2 



1 1 

dz 
1 



dx 



(l-x) N 



J2e( s (z-~z)) 

s=± 



dtt N ~ 



[F(mrz-,P,p)] N 

In addition we make t— dependence of F(Q(t), z;p, P) explicit 



where 



F(fl(t) , 5; p, P) = t + F(fi(0), 5; p, P) . 

1 — x 



1 — x x 
s T/ x b + a b — a _ 

P(t) = J(0J,z)t + — ; —2, 



J(o», 2) = 



2 2 
6 + a 6 — a 



(A25) 



(A26) 



(A27) 



(A28) 



Integrate over the variable t we get 
dtt N ~ 2 



1 



t n 



F(Sl(t), z; p, P)" (N - 1) F(SI(0), z; p, P) [F(Q(t), ~z; p, P)] 



N-l 



(A29) 
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and hence 



2(4tt) 



/ dz dx~S^ 

J-i Jo t± (1 - 



6[s(z - ~z)]T t 



N-l 
s 



x) N F{tt{0), z; p, P) [F(A(T s ),z; p, P)} 



N-l 



(A30) 



In order to separate F's in the denominator we use the following identity 

1 



F(O(0),z;p,P) [F(n(t),z;p,P)} 
1 — x 
J(u, z)T s 



N-l 



F(n(0) } z;p,P) F(A(T s ),z;p,P) 



1 



[F(A(T s ),z;p,P)} 



N-2 



(A31) 



Note, for a given a given N one can repeat N-l times this algebra until the power of the 
last factor vanishes, which is just the reason of introducing the operation (IA31I) . After this 
operation the momentum dependence of the denominators at each term turns to be formally 
the same as in the desired IR. Although it is possible to derive the appropriate formula for 
any N, it would lead to an unpractical expressions (probably it cannot be written in a closed 
form) and we rather choose one concrete value of N. Motivated by the success of the scalar 
model studies, we take N = 2 from now. Explicitly written we have for I 



1 



2(4yrr J-i 



/ dz I dx Y"* 

J-i Jo £± 



0[s(z ~ z)] 
J(u, z){l — x) 



1 



F(n(0),z;p,P) F(A(T 8 ),z;p,P) 



(A32) 



Integrating per-partes in variable x we get for f)A32|) 



-1 



2(4tt) 



dz dx 



9[s(z - z)] ln(l - x) 



s=± 



J(u, z) 



dQ,(0) 
dx 



dA(T a ) 
dx 



F 2 (fi(0),5;p,P) F\A{T s )rz;p,P)_ 

(A33) 

Implementing the following identity into the integrand of ( 1A33I) 



5{u - Q(t)) 



(A34) 



into the integrand of (1A33I) and changing the order of the integration and performing inte- 
gration in the variable x, we arrive to the desired result for the second term in the second 
line of the Rel. Ml : 



oo rl 

duj I dz 

J -i [F(uj,z;p,P)) 



N ' 



X3{U,Z) 



Pv(a)p v (b) ^ / ln(l -^(0)) 



s 2(4tt) 2 



3=± 



J(u) } z) 



sgn 



53 d [s(z - S)] Hl T/ X ^ s)) 9[D(0)]0[x 3 (T s )]e[l - x^sga {E[ Xj (T.)]} 



s=± 



J(oo, z) 



dxj(0) 

(A35) 
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where we have included omitted prefactor in (IA9I) and Xj are the roots of the delta function 
argument in (1A34|) 



X±(T) = 

D(T) = B(T) 2 - AAc , B(T) = R(T) - u - c 

dtt(t) E(x) „, . c . , . 

— ^ = , E(x) =u-S--. (A36) 

da; 1 — x x A 

Introducing a bit compact notation for the sum over an arbitrary function U with pa- 
rameter T 

E U{T) = U(0) - 9(z - z)U(T + ) - 6(~z - z)U(T_) (A37) 

T 

we can rewrite the above relation (IA35I) like 

xs&z) = I E E H )~ { *f )] e[D{T)]o[*s(F)W - xsW^Mn- 

(A38) 

The first term in the second line of Rel. flA9j) should be derived can treated in very similar 
fashion as the previous case, note only a different power of F in the denominator. Doing 
this explicitly and adding a correct prefactor the appropriate expression reads: 



2 



/ p2 a _ b \ 

p v (a)p v (b) ( u - — + — ^— z I - p s (a)p s (b) 



(A39) 

Assuming the validity of Theorem o the Uniqueness then the momentum BSE is converted 
to the equation for the weight function. It reads 



p [2] (Co,z)= [ dco [ dzV [2] (Co,z;uo,z)p [2] (uj,z) , (A40) 



• i 

10 J-l 

where the kernel is simply given by the sum of contributions derived above 

pco poo poo 4 

V*- 2 ' (Co, z; uo, z) = — 3g 2 / del da I dbpJc)^ Xi(^iZ\uo, z). (A41) 

Jo Jo Jo ~{ 

Heavy quark approximation -unequal mass case 

When the quark is sufficiently heavy (say m q (2GeV) » Aqcd) then the approximation 
of quark mass function by a constant appears to be adequate. Neglecting the selfenergy 
correction is equivalent to use a free heavy quark propagator, this corresponds with the free 
particle spectral functions: 

M lPv (a) = p s (a) = M l 8(a - M 2 ) , 

M 2 p v (b) = Ps (b) = M 2 5(b - M 2 ) , (A42) 
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where the variable a, b distinguish the type of quark from which the bound state is composed. 
For completeness we write down the kernel explicitly here for this case. 

3q 2 c°° 

V [2] (a), z; w, z) = — ^- / dcp g (c) ( X i + X2 + X3 + Xa) I (A43) 



2(4tt) 2 Jo 



z9(z -z)^ y a ±9(D a )e(y a± )0(l - y a± ) 



Xl T+ l + z 5 (i-y a ±)VD~ a 

zO(z - z) y b ±6(D b )6(y b ±)6(l - y b± ) 
X2 = T— — ^ 



± (l-y b ±)VD b 
X 4 = f oj - — + 2 2 z - M X M 2 \ 

?s g^ ( _^ - ln(1 _ , ATMT)] } 

(A44) 



and where the arguments a, b in the functions x, J, D and must be replaced by the quark 
masses M ± , M 2 . 

Equal mass case 

In the case of quarkonia the kernel becomes more symmetric with respect to the variable z 
and the formula for the kernel further simplify. The function R depends on z l s only through 
the variable T such that 

J(u,z) -> J = u - M 2 ; R(T) -> R(T) = JT + M 2 (A45) 

where M is a common mass M = Mi = M 2 . The roots become 

y a ± -> Xj=±(7+); y 6 ± -> Xj=±(T_) , (A46) 

in this case and the kernel can be written in more compact form: 

VW(cu,z;u; : z) = J™ dcp g (c)K^(Cu, z;u, z, c) (A47) 

K^(uj, z; u, z, c) = 

E E - 5)]^(Tj]^(T s )]tf[l - Xj (T s )}- ~ szT ^ T s) 

s=±l=± (1 + sz)y/D(T 8 )(l - Xj(T s )) 

r j=± 

— fa^ - - -ln(l -^.(T))sgni^(T)]^- }> (A l,s) 



|^(T)]|(1-^(T)) v //o , J2J2 
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One gluon exchange approximation 

In order to consider one gluon approximation we should take the massless limit c — > 0. 
In addition we also restrict ourselves to the equal mass case. One can easily recognize that 
the root (x + = 0) is trivial and associated contribution in flA480 vanishes. 

For the purpose of brevity we label 

x(T)=X-(T) = *: Ri p . (A49) 
Taking into account the relations 

Vd/x = a,e = a, 

and doing a little algebra we get for the kernel: 

rr\2] , v ~3g 2 f P 2 n , 2X1 (m 2 -S\ 

V OGE (u;,z;u,z) = — j— 0(u, - m ) In [—j) 
+ Y,6[s{z-z)]6(-B a )6(A + B a 

s=± 



A + B s V A J AJ 2 



For completeness we review the complete list of the functions at this place: 

p2 

A = Cu-S, J = u-M 2 ,S=(l-z 2 ) — 

B s = (uu- M 2 )T S + M 2 -u, T± = — ^. (A51) 

-L _L_ 2 

Case N=l for 

Repeating the derivation for the parameter iV = 1 we should obtain the following homo- 
geneous equation 

poo pi 

p [1] (u,z) = J^ duo J ^dzV [1] ( y u,z;u,z)p [1 \u,z) (A52) 
where the kernel is given by the following expression 



V [1] (Co,z;u,z) 



s 2(4tt) 2 ' 



Xi = -p v (a)p v {b)——0(z -z)Xj 7 =^= J -. 

1 — z - =± V L>b 

g[D(0)]6)[x 3 -(0)]g[l-x 3 -(0)] 



j=± ~ L*J < 



|£[x,(T)]|J(^) 
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where we have used the same notations and conventions as in the previous Section. The 
appropriate derivation is very straightforward and exactly repeats the steps made for the 
case with N = 2, hence we make only few comments here. The relations for xi,2 are in fact 
derived in the previous section (since for general N) . The relation for \i is adopted from the 



paper lljand remaining function X3 follows from the conversion of the term with F° = 1, 
i.e. 

A 1 (A54) 



(27r) 4 £>i£> 2 iV 

which is presented in the momentum space BSE in that case. The appropriate derivation 
is in fact more easy then in the case N = 2 and the appropriate result represents the basic 
Perturbation Theory Integral Representation for a scalar triangle with one leg momentum 
constrained such that P 2 < (M x + M 2 ) 2 . 

APPENDIX B: NUMERICAL PROCEDURE 

In this section we describe the numerical procedure actually used for obtaining the spectra 
of the bound states. Because of the structure of the integral equations to be solved, the 
treatment is similar the procedure used in the papers [Io|, however we have used some 



tricky improvements modification here which enforce the numeric stability when the mass 
parameter c in gluon IR (or simple vector boson mass m c ) is to small compared to the 
constituent mass. We describe the technical difference in that case. As a first we describe 
the numerical treatment for the case c ~ m. 

The equation (IA48I) represents the homogeneous linear integral equation which solution 
needs to be properly normalized. Adopting the following auxiliary normalization condition 
for this purpose: 



1 Jo 



1= dz du V/ ' . (Bl) 



J 2 



Then we can find that the equal constituent mass BSE can be transformed into the inho- 
mogeneous integral equation (1B2j) 



p [2] (Cu, ~z) = Kf (u, ~z) + duj J 1 dzK [ H ] (Co, z; u, z)p [2] (co, z) (B2) 



with the following kernels: 

K?(u,, ~z) = ^0[£>(O)M%(O)]0[1 - ,x,(0)] ln(l - x,(0))sgnE[x,(0)], (B3) 

K$(Co,z;u,z,c) = g e[D(0)]9[ Xj (0)W -^(0)]ln(l - Xj ( )) ^^M - 1 

In 2 

+ tJ~T 2 E E 0[s{z - ~z)}e[D(T s )]e\x 3 (T s )]6\l - Xj (T s )} 

z l 47r J s=±j=± 

szT s Xj(T s ) 



l + sz)JD(T s )(l- Xj (T s )) 
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(i-S)t. +! »< L ^^ £. (B4) 



1^(^)11(1 -x^T,)) J yJ 2J 

The meaning of the all functions were established in the previous appendix, for convenience 
of the reader we review them also here: 



-B(T) ± JD(T) 
x±(T) = J 2A y ; D(T) = B(Tf-AAm c , 

A = Cu-S; B(T) = R(T) -Co~m Cl 

E(x) = uj-^-S; R(T) = JT + m 2 ; J = u - m 2 , 
x 2 

J = tu-m 2 ; S=(l-~z 2 )^-; T± = . (B5) 

The equation ( 1B2I) has been actually used for the numerical solution when m c /m q ~ 1. 

The kernel t]B2[) is free of the running singularities because of the presence of theta 
functions, however in the case of massless gluon kernel we would have the singularity just 
on the boundary. There J — > as u approaches the quarks threshold. We have found that 
this instability is avoided if we generate the inhomogeneous term in a following manner: We 
formally add the zero in the following form 

(f(Co,z)-mz))^^ (B6) 

ur 

to the right hand side of the original homogeneous equation Bethe-Salpeter equation for 
weight function p. Then solving the equation 



p' 2 '(cc>, z) = f(u, z) + J duo J dz 
within the following normalization condition 



z;u,z 



p®(u,,z) (B7) 



1 dz r d j w ^\ ( B8) 

-1 JO IjJ 1 

is equivalent to the solution of the original BSE. As a suitable function we choose 

f(u,z) = l (B9) 

We have found that this method is applicable for any positive value of m 2 , but is slowly 
convergent when is used for previously discussed case m c ~ m q . Note here, that up to the 
small numerical error the equation (1B2I) offers the same spectra for such a's and m! g s where 
the both equations ( 1B2I) . (1B7I) are numerically stable. 

The equation (1B7I) has been solved by the method of the iteration. If the iterations 
failed- measure being both the difference of the rhs. and lhs. of the integral equation and 
deviation of the auxiliary normalization integral from a predefined value- we were changing 
the coupling constant (in the treatment with P 2 fixed, otherwise procedure is opposite) until 
the solution was found. For numerical solution we discretize integration variables to and z 
using Gauss-Legendre quadratures (with suitable mapping from < — 1, +1 >^< LU min , oo > 
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for u ). Equations (1B7I) . (1B2I) is solved on the grid of N = N z * points which are spread 
on the rectangle (—1, +1) * (u) m i n , oo). The value ui m i n is given by the support of the spectral 
function. In the all cases we take N u = 2N Z . Examples of a numerical convergence for 
some cases of bound states are presented in Table III. As we can see, there is a rather 
weak dependence of the eigenvalue a on the number of mash points N^. The last value is 
calculated from the weighted average (WA) with N u be the appropriate weight. 





32 


40 


64 


80 


WA 


r? = 0.95; m g /M = 10~ 3 


0.6611 


0.6690 


0.6697 


0.6734 


0.669 


•q = 0.95; m g /M = 0.5 


1.037 


1.0259 


1.0210 


1.0229 


1.029 


r] = 0.99; m g /M = 0.5 


0.818 


0.8127 


0.8158 


0.8155 


0.816 



TABLE 3. The coupling a s = g 2 /(4ir) for the ladder BSE with fixed ratio m c /M as a 

function of the number of mesh-points. 
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